﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Text;
using System.Threading.Tasks;

namespace LightCAD.MathLib
{
    internal class EllipseAndEllipse
    {
        public static List<double> SolveQuadraticEquation(double a,double b,double c)
        {
            List<double> root = new List<double>();
            b = b / a;
            c = c / a;
            double D = b * b - 4 * c;
            if (Math.Abs(D) < 1e-7)
            {
                root.Add((-0.5 * b));
            }
            if (D < 0)
            {

            }
            double sqrtD = Math.Sqrt(D);
            root.Add((0.5 * (-b - sqrtD)));
            root.Add((0.5 * (-b + sqrtD)));
            return root;
        }

        public static List<double> SolveQuarticEquation(double a, double b, double c,double d)
        {
            double p = b - (a * a * 3 / 8);
            double q =  a * a * a * 1 / 8 - 0.5 * a * b + c;
            double r = d - 0.25 * a * c +   a * a * b * 1 / 16 -  a * a * a * a * 3 / 256;

            List<double> root = new List<double>();

            if (Math.Abs(q) < 1e-7)
            {
                double D = p * p - 4 * r;
                if (Math.Abs(D) < 1e-5)
                {
                    double m = -0.5 * p;
                    if (m >= 0)
                    {
                        root.Add(Math.Sqrt(m));
                    }
                    if (m > 0)
                    {
                        root.Add(-Math.Sqrt(m));
                    }
                }
                else if (D > 0)
                {
                    double sqrtD = Math.Sqrt(D);
                    double m1 = (-p - sqrtD) * 0.5;
                    double m2 = (-p + sqrtD) * 0.5;
                    if (m1 >= 0)
                    {
                        root.Add(Math.Sqrt(m1));
                    }
                    if (m1 > 0)
                    {
                        root.Add(-Math.Sqrt(m1));
                    }
                    if (m2 >= 0)
                    {
                        root.Add(Math.Sqrt(m2));
                    }
                    if (m2 > 0)
                    {
                        root.Add(-Math.Sqrt(m2));
                    }
                }
            }
            else
            {
                Vector2 t = GetOneCubicEquationRoot(2 * p, p * p - 4 * r, -q * q);
                Vector2 z = sqrt(t);
                Vector2 u = Mul(new Vector2(p, 0).Add(t),0.5);
                Vector2 v=Mul(inverse(z), q * 0.5);
                Vector2 u2 = new Vector2(u.X, u.Y);
                Vector2 u3 = new Vector2(u.X, u.Y);
                List<Vector2> x12 = SolveComplexQuadraticEquation(z, u2.Sub(v));
                List<Vector2> x34 = SolveComplexQuadraticEquation(Mul(z,-1), u3.Add(v));

                if (Math.Abs(x12[0].Y) < 1e-7 )
                {
                    root.Add(x12[0].X);
                }
                if (Math.Abs(x12[1].Y) < 1e-7)
                {
                    root.Add(x12[1].X);
                }
                if (Math.Abs(x34[0].Y) < 1e-7)
                {
                    root.Add(x34[0].X);
                }
                if (Math.Abs(x34[1].Y) < 1e-7)
                {
                    root.Add(x34[1].X);
                }
                   
            }
            for (var i = 0; i < root.Count; i++)
            {
                root[i] -= 0.25 * a;
            }
            return root;
           
        }
        public static Vector2 GetOneCubicEquationRoot(double a,double b,double c)
        {
            double p = b - a * a * 0.33333333333;
            double q =   a * a * a * 2 / 27 - a * b * 0.33333333333 + c;
            Vector2 l =  new Vector2(0.03703703703 * p * p * p + 0.25 * q * q, 0);
            Vector2 sqrtQ = sqrt(l);
            Vector2 A = Qbrt(new Vector2(-0.5 * q).Sub(sqrtQ));
            if (A.X == 0 && A.Y == 0)
            {
                A = Qbrt(new Vector2(-0.5 * q).Add(sqrtQ));
            }
            Vector2 B = Mul(inverse(A), -p * 0.33333333333);
            return A.Add(B).Sub(new Vector2(a * 0.33333333333));
        }
        public static List<Vector2> SolveComplexQuadraticEquation(Vector2 b, Vector2 c)
        {
            Vector2 sqrtD1 =new Vector2((b.X * b.X - b.Y * b.Y) - c.X * 4,
                (b.X * b.Y + b.X * b.Y) - c.Y * 4);
            Vector2 sqrtD = new Vector2();
            sqrtD= sqrt(sqrtD1);
            List<Vector2> list = new List<Vector2>();
            Vector2 point1 = new Vector2((b.X + sqrtD.X) * (-0.5), (b.Y + sqrtD.Y) * (-0.5));
            Vector2 point2 = new Vector2((b.X - sqrtD.X) * (-0.5), (b.Y - sqrtD.Y) * (-0.5));
            list.Add(point1);
            list.Add(point2);
            return list;


        }
        public static Vector2 sqrt (Vector2 b)
        {
            double l = Math.Sqrt(b.X * b.X + b.Y * b.Y);
            double sgn = 0;
            if (b.Y > 0)
            {
                sgn = 1;
            }
            if(b.Y == 0)
            {
                sgn = 0;
            }
            if (b.Y < 0)
            {
                sgn = -1;
            }
            if (sgn == 0)
            {
                sgn = 1;
            }
            return new Vector2(Math.Sqrt(0.5*(l+b.X)), sgn*Math.Sqrt(0.5*(l-b.X)));
        }

        public static Vector2 Qbrt(Vector2 complex)
        {
            double angle = Math.Atan2(complex.Y, complex.X) * 0.33333333333;
            double l = Math.Pow(complex.X * complex.X + complex.Y * complex.Y, 0.16666666666);
            return new Vector2(l * Math.Cos(angle), l * Math.Sin(angle));
        }

        public static Vector2 inverse(Vector2 a)
        {
            double  l = 1.0 / (a.X * a.X + a.Y * a.Y);
            return new Vector2(a.X * l, -a.Y * l);
        }

        public static Vector2 Mul(Vector2 a,double factor)
        {
            return new Vector2(a.X * factor, a.Y * factor);
        }
    }
}
